#!/usr/bin/env Rscript
##################################
############ Preamble ############
##################################

# set language to English
Sys.setenv(LANG = "en")

# clean up
rm(list = ls())

# Set working directory: Please set your own
if (Sys.getenv("RSTUDIO") == "1") setwd("~/Dropbox/cues_bjps_replication")

# Load necessary packages
library(Rmisc)
library(tidyverse)
library(readxl)
library(stargazer)
library(kableExtra)
library(sandwich)
library(broom)
library(dotwhisker)
library(jtools)
library(scales)
library(xtable)

# read data
cols <- as.character(read_excel("data/i360_dataset_Normwandel_220201.xlsx", n_max = 1, col_names = FALSE))
raw_data <- read_excel("data/i360_dataset_Normwandel_220201.xlsx", skip = 2, col_names = cols)


# rename variables
raw_data <- raw_data %>%
  dplyr::rename(
    disposition_code = dispcode,
    treatment = c_0031,
    cover_page = v_27,
    age = v_28,
    age_range = v_29,
    gender = v_30,
    state = v_31,
    born_in_germany = v_39,
    sign_petition_self = v_47,
    sensitive_item_agree = dupl1_v_44,
    left_right = dupl1_v_41,
    petition_appropriate_self = v_48,
    sign_petition_others = v_49,
    petition_appropriate_others = v_50,
    delete_tweet = v_51,
    voted_2021 = v_53,
    party_voted_2021 = v_54,
    party_voted_2021_other_which = v_55,
    voted_2017 = v_56,
    party_voted_2017 = v_57,
    party_voted_2017_other_which = v_58,
    household_income = v_59,
    education = v_60,
    mother_born_where = v_61,
    mother_born_other_where = v_62,
    father_born_where = v_63,
    father_born_other_where = v_64
  )

# Fix treatment arm labels
raw_data <- raw_data %>%
  mutate(treatment = replace(treatment, which(treatment == 1), "Mainstream Approve")) %>%
  mutate(treatment = replace(treatment, which(treatment == 2), "Mainstream Approve")) %>%
  mutate(treatment = replace(treatment, which(treatment == 3), "RRP Approve")) %>%
  mutate(treatment = replace(treatment, which(treatment == 4), "RRP Approve")) %>%
  mutate(treatment = replace(treatment, which(treatment == 5), "Mainstream and RRP Approve")) %>%
  mutate(treatment = replace(treatment, which(treatment == 6), "Mainstream and RRP Approve")) %>%
  mutate(treatment = replace(treatment, which(treatment == 7), "Mainstream Disapprove and RRP Approve")) %>%
  mutate(treatment = replace(treatment, which(treatment == 8), "Mainstream Disapprove and RRP Approve")) %>%
  mutate(treatment = replace(treatment, which(treatment == 9), "Control")) %>%
  mutate(treatment = replace(treatment, which(treatment == 10), "Control"))

# recode state names
raw_data <- raw_data %>%
  mutate(state = replace(state, which(state == 1), "Baden-Wuerttemberg")) %>%
  mutate(state = replace(state, which(state == 2), "Bayern")) %>%
  mutate(state = replace(state, which(state == 3), "Berlin")) %>%
  mutate(state = replace(state, which(state == 4), "Brandenburg")) %>%
  mutate(state = replace(state, which(state == 5), "Bremen")) %>%
  mutate(state = replace(state, which(state == 6), "Hamburg")) %>%
  mutate(state = replace(state, which(state == 7), "Hessen")) %>%
  mutate(state = replace(state, which(state == 8), "Mecklenburg-Vorpommern")) %>%
  mutate(state = replace(state, which(state == 9), "Niedersachsen")) %>%
  mutate(state = replace(state, which(state == 10), "Nordrhein-Westfalen")) %>%
  mutate(state = replace(state, which(state == 11), "Rheinland-Pfalz")) %>%
  mutate(state = replace(state, which(state == 12), "Saarland")) %>%
  mutate(state = replace(state, which(state == 13), "Sachsen")) %>%
  mutate(state = replace(state, which(state == 14), "Sachsen-Anhalt")) %>%
  mutate(state = replace(state, which(state == 15), "Schleswig-Holstein")) %>%
  mutate(state = replace(state, which(state == 16), "Thueringen"))

# recode age_ranges
raw_data <- raw_data %>%
  mutate(age_range = replace(age_range, which(age_range == 1), "18 - 29 years")) %>%
  mutate(age_range = replace(age_range, which(age_range == 2), "30 - 39 years")) %>%
  mutate(age_range = replace(age_range, which(age_range == 3), "40 - 49 years")) %>%
  mutate(age_range = replace(age_range, which(age_range == 4), "50 - 59 years")) %>%
  mutate(age_range = replace(age_range, which(age_range == 5), "60 - 90 years"))

# recode gender
raw_data <- raw_data %>%
  mutate(gender = replace(gender, which(gender == 1), "Male")) %>%
  mutate(gender = replace(gender, which(gender == 2), "Female"))

# recode sensitive item - change sensitive item response "disagree" to 0
raw_data <- raw_data %>%
  mutate(sensitive_item_agree = replace(sensitive_item_agree, which(sensitive_item_agree == 2), 0))

# recode coordination game - change "not willing to sign" from 2 to 0
raw_data <- raw_data %>%
  mutate(sign_petition_self = replace(sign_petition_self, which(sign_petition_self == 2), 0))

# recode punishment exp - change "do not delete tweet" from 2 to 0
raw_data <- raw_data %>%
  mutate(delete_tweet = replace(delete_tweet, which(delete_tweet == 2), 0))

# recode left_right 100 is 0, 99 is NA
raw_data <- raw_data %>%
  mutate(left_right = replace(left_right, which(left_right == 100), 0)) %>%
  mutate(left_right = replace(left_right, which(left_right == 99), NA))

# recode income: recode "don't know" responses to NA
raw_data <- raw_data %>%
  mutate(household_income = replace(household_income, which(household_income == 5), NA))

# reverse scale income because it is inverted (currently 1 is highest, 4 is lowest)
raw_data <- raw_data %>% mutate(household_income = 5 - household_income)

# recode no school leaving certificate to 1 as it is currently 9 and the highest in the scale, and move all other values in scale up by 1
raw_data <- raw_data %>%
  mutate(education = replace(education, which(education == 9), 0))

raw_data <- raw_data %>% mutate(education = 1 + education)

# recode party voted
raw_data <- raw_data %>%
  mutate(party_voted_2021 = replace(party_voted_2021, which(party_voted_2021 == 1), "CDU/CSU")) %>%
  mutate(party_voted_2021 = replace(party_voted_2021, which(party_voted_2021 == 2), "SPD")) %>%
  mutate(party_voted_2021 = replace(party_voted_2021, which(party_voted_2021 == 3), "AfD")) %>%
  mutate(party_voted_2021 = replace(party_voted_2021, which(party_voted_2021 == 4), "Green")) %>%
  mutate(party_voted_2021 = replace(party_voted_2021, which(party_voted_2021 == 5), "FDP")) %>%
  mutate(party_voted_2021 = replace(party_voted_2021, which(party_voted_2021 == 6), "Die Linke"))

# remove incompletes or inattentives, retain only those who completed the survey or took a break but finished the survey
analysis_data <- raw_data %>%
  filter(disposition_code == "31" | disposition_code == "32")

# clean up sensitive item - drop -77s, merge both columns
analysis_data <- analysis_data %>%
  mutate(sensitive_item_agree = replace(sensitive_item_agree, which(sensitive_item_agree < 0), NA))

# Create dummy variable
analysis_data$MRPapprovedummy <- ifelse(analysis_data$treatment == "Mainstream Approve", 1, 0)
analysis_data$RRPapprovedummy <- ifelse(analysis_data$treatment == "RRP Approve", 1, 0)
analysis_data$MRPapproveRRPapprovedummy <- ifelse(analysis_data$treatment == "Mainstream and RRP Approve", 1, 0)
analysis_data$MRPdisapproveRRPapprovedummy <- ifelse(analysis_data$treatment == "Mainstream Disapprove and RRP Approve", 1, 0)

# right-wing only (including AfD)
right_only <- analysis_data %>%
  filter(party_voted_2021 == "CDU/CSU" | party_voted_2021 == "FDP" | party_voted_2021 == "AfD")

# left-wing only
left_only <- analysis_data %>%
  filter(party_voted_2021 == "SPD" | party_voted_2021 == "Green" | party_voted_2021 == "Die Linke")

# standardize by subtracting mean and dividing by SD
analysis_data[, c(16, 17, 18, 19, 20, 21)] <- scale(analysis_data[, c(16, 17, 18, 19, 20, 21)])
left_only[, c(16, 17, 18, 19, 20, 21)] <- scale(left_only[, c(16, 17, 18, 19, 20, 21)])
right_only[, c(16, 17, 18, 19, 20, 21)] <- scale(right_only[, c(16, 17, 18, 19, 20, 21)])

# Check number of observations in each subsample
table(left_only$treatment)
table(right_only$treatment)

# regression for Agreement with Sensitive Item across different samples
sensitive_item_nocontrols_fullsample <- lm(sensitive_item_agree ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = analysis_data)
sensitive_item_nocontrols_leftonly <- lm(sensitive_item_agree ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = left_only)
sensitive_item_nocontrols_rightonly <- lm(sensitive_item_agree ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = right_only)

# robust SEs - "HC1" for STATA
robust_sensitive_item_nocontrols_fullsample <- summ(sensitive_item_nocontrols_fullsample, robust = "HC1")
robust_sensitive_item_nocontrols_leftonly <- summ(sensitive_item_nocontrols_leftonly, robust = "HC1")
robust_sensitive_item_nocontrols_rightonly <- summ(sensitive_item_nocontrols_rightonly, robust = "HC1")

# convert regression outcomes to df for coefficient plot
robust_sensitive_item_nocontrols_fullsample_df <- broom::tidy(robust_sensitive_item_nocontrols_fullsample) %>%
  mutate(sample = "Full Sample") %>%
  mutate(measure = "Agreement with Sensitive Item") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sensitive_item_nocontrols_fullsample_df)[1] <- "model"
names(robust_sensitive_item_nocontrols_fullsample_df)[7] <- "term"

robust_sensitive_item_nocontrols_leftonly_df <- broom::tidy(robust_sensitive_item_nocontrols_leftonly) %>%
  mutate(sample = "Left Wing Only") %>%
  mutate(measure = "Agreement with Sensitive Item") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sensitive_item_nocontrols_leftonly_df)[1] <- "model"
names(robust_sensitive_item_nocontrols_leftonly_df)[7] <- "term"


robust_sensitive_item_nocontrols_rightonly_df <- broom::tidy(robust_sensitive_item_nocontrols_rightonly) %>%
  mutate(sample = "Right Wing Only") %>%
  mutate(measure = "Agreement with Sensitive Item") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sensitive_item_nocontrols_rightonly_df)[1] <- "model"
names(robust_sensitive_item_nocontrols_rightonly_df)[7] <- "term"

# Willingness to Sign Petition

# regression for Willingness to Sign Petition across different samples
sign_petition_self_nocontrols_fullsample <- lm(sign_petition_self ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = analysis_data)
sign_petition_self_nocontrols_leftonly <- lm(sign_petition_self ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = left_only)
sign_petition_self_nocontrols_rightonly <- lm(sign_petition_self ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = right_only)

# robust SEs
robust_sign_petition_self_nocontrols_fullsample <- summ(sign_petition_self_nocontrols_fullsample, robust = "HC1")
robust_sign_petition_self_nocontrols_leftonly <- summ(sign_petition_self_nocontrols_leftonly, robust = "HC1")
robust_sign_petition_self_nocontrols_rightonly <- summ(sign_petition_self_nocontrols_rightonly, robust = "HC1")

# convert regression outcomes to df for coefficient plot
robust_sign_petition_self_nocontrols_fullsample_df <- broom::tidy(robust_sign_petition_self_nocontrols_fullsample) %>%
  mutate(sample = "Full Sample") %>%
  mutate(measure = "Personal Willingness to Sign Petition") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sign_petition_self_nocontrols_fullsample_df)[1] <- "model"
names(robust_sign_petition_self_nocontrols_fullsample_df)[7] <- "term"

robust_sign_petition_self_nocontrols_leftonly_df <- broom::tidy(robust_sign_petition_self_nocontrols_leftonly) %>%
  mutate(sample = "Left Wing Only") %>%
  mutate(measure = "Personal Willingness to Sign Petition") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sign_petition_self_nocontrols_leftonly_df)[1] <- "model"
names(robust_sign_petition_self_nocontrols_leftonly_df)[7] <- "term"


robust_sign_petition_self_nocontrols_rightonly_df <- broom::tidy(robust_sign_petition_self_nocontrols_rightonly) %>%
  mutate(sample = "Right Wing Only") %>%
  mutate(measure = "Personal Willingness to Sign Petition") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sign_petition_self_nocontrols_rightonly_df)[1] <- "model"
names(robust_sign_petition_self_nocontrols_rightonly_df)[7] <- "term"

# Personal Views about Appropriateness of Signing Petition
# regression for Personal Views about Appropriateness of Signing Petition across different samples
petition_appropriate_self_nocontrols_fullsample <- lm(petition_appropriate_self ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = analysis_data)
petition_appropriate_self_nocontrols_leftonly <- lm(petition_appropriate_self ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = left_only)
petition_appropriate_self_nocontrols_rightonly <- lm(petition_appropriate_self ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = right_only)

# robust SEs
robust_petition_appropriate_self_nocontrols_fullsample <- summ(petition_appropriate_self_nocontrols_fullsample, robust = "HC1")
robust_petition_appropriate_self_nocontrols_leftonly <- summ(petition_appropriate_self_nocontrols_leftonly, robust = "HC1")
robust_petition_appropriate_self_nocontrols_rightonly <- summ(petition_appropriate_self_nocontrols_rightonly, robust = "HC1")

# convert regression outcomes to df for coefficient plot
robust_petition_appropriate_self_nocontrols_fullsample_df <- broom::tidy(robust_petition_appropriate_self_nocontrols_fullsample) %>%
  mutate(sample = "Full Sample") %>%
  mutate(measure = "Personal Views of Appropriateness of Signing") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_petition_appropriate_self_nocontrols_fullsample_df)[1] <- "model"
names(robust_petition_appropriate_self_nocontrols_fullsample_df)[7] <- "term"

robust_petition_appropriate_self_nocontrols_leftonly_df <- broom::tidy(robust_petition_appropriate_self_nocontrols_leftonly) %>%
  mutate(sample = "Left Wing Only") %>%
  mutate(measure = "Personal Views of Appropriateness of Signing") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_petition_appropriate_self_nocontrols_leftonly_df)[1] <- "model"
names(robust_petition_appropriate_self_nocontrols_leftonly_df)[7] <- "term"


robust_petition_appropriate_self_nocontrols_rightonly_df <- broom::tidy(robust_petition_appropriate_self_nocontrols_rightonly) %>%
  mutate(sample = "Right Wing Only") %>%
  mutate(measure = "Personal Views of Appropriateness of Signing") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_petition_appropriate_self_nocontrols_rightonly_df)[1] <- "model"
names(robust_petition_appropriate_self_nocontrols_rightonly_df)[7] <- "term"


# Empirical Expectations
# regression for Empirical Expectations across different samples
sign_petition_others_nocontrols_fullsample <- lm(sign_petition_others ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = analysis_data)
sign_petition_others_nocontrols_leftonly <- lm(sign_petition_others ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = left_only)
sign_petition_others_nocontrols_rightonly <- lm(sign_petition_others ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = right_only)

# robust SEs
robust_sign_petition_others_nocontrols_fullsample <- summ(sign_petition_others_nocontrols_fullsample, robust = "HC1")
robust_sign_petition_others_nocontrols_leftonly <- summ(sign_petition_others_nocontrols_leftonly, robust = "HC1")
robust_sign_petition_others_nocontrols_rightonly <- summ(sign_petition_others_nocontrols_rightonly, robust = "HC1")

# convert regression outcomes to df for coefficient plot
robust_sign_petition_others_nocontrols_fullsample_df <- broom::tidy(robust_sign_petition_others_nocontrols_fullsample) %>%
  mutate(sample = "Full Sample") %>%
  mutate(measure = "Empirical Expectations") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sign_petition_others_nocontrols_fullsample_df)[1] <- "model"
names(robust_sign_petition_others_nocontrols_fullsample_df)[7] <- "term"


robust_sign_petition_others_nocontrols_leftonly_df <- broom::tidy(robust_sign_petition_others_nocontrols_leftonly) %>%
  mutate(sample = "Left Wing Only") %>%
  mutate(measure = "Empirical Expectations") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sign_petition_others_nocontrols_leftonly_df)[1] <- "model"
names(robust_sign_petition_others_nocontrols_leftonly_df)[7] <- "term"


robust_sign_petition_others_nocontrols_rightonly_df <- broom::tidy(robust_sign_petition_others_nocontrols_rightonly) %>%
  mutate(sample = "Right Wing Only") %>%
  mutate(measure = "Empirical Expectations") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_sign_petition_others_nocontrols_rightonly_df)[1] <- "model"
names(robust_sign_petition_others_nocontrols_rightonly_df)[7] <- "term"

# Normative Expectations
# regression for Normative Expectations across different samples
petition_appropriate_others_nocontrols_fullsample <- lm(petition_appropriate_others ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = analysis_data)
petition_appropriate_others_nocontrols_leftonly <- lm(petition_appropriate_others ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = left_only)
petition_appropriate_others_nocontrols_rightonly <- lm(petition_appropriate_others ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = right_only)

# robust SEs
robust_petition_appropriate_others_nocontrols_fullsample <- summ(petition_appropriate_others_nocontrols_fullsample, robust = "HC1")
robust_petition_appropriate_others_nocontrols_leftonly <- summ(petition_appropriate_others_nocontrols_leftonly, robust = "HC1")
robust_petition_appropriate_others_nocontrols_rightonly <- summ(petition_appropriate_others_nocontrols_rightonly, robust = "HC1")

# convert regression outcomes to df for coefficient plot
robust_petition_appropriate_others_nocontrols_fullsample_df <- broom::tidy(robust_petition_appropriate_others_nocontrols_fullsample) %>%
  mutate(sample = "Full Sample") %>%
  mutate(measure = "Normative Expectations") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_petition_appropriate_others_nocontrols_fullsample_df)[1] <- "model"
names(robust_petition_appropriate_others_nocontrols_fullsample_df)[7] <- "term"

robust_petition_appropriate_others_nocontrols_leftonly_df <- broom::tidy(robust_petition_appropriate_others_nocontrols_leftonly) %>%
  mutate(sample = "Left Wing Only") %>%
  mutate(measure = "Normative Expectations") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_petition_appropriate_others_nocontrols_leftonly_df)[1] <- "model"
names(robust_petition_appropriate_others_nocontrols_leftonly_df)[7] <- "term"

robust_petition_appropriate_others_nocontrols_rightonly_df <- broom::tidy(robust_petition_appropriate_others_nocontrols_rightonly) %>%
  mutate(sample = "Right Wing Only") %>%
  mutate(measure = "Normative Expectations") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_petition_appropriate_others_nocontrols_rightonly_df)[1] <- "model"
names(robust_petition_appropriate_others_nocontrols_rightonly_df)[7] <- "term"

# Sanctioning
# regression for Sanctioning across different samples
delete_tweet_nocontrols_fullsample <- lm(delete_tweet ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = analysis_data)
delete_tweet_nocontrols_leftonly <- lm(delete_tweet ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = left_only)
delete_tweet_nocontrols_rightonly <- lm(delete_tweet ~ MRPapprovedummy + RRPapprovedummy + MRPdisapproveRRPapprovedummy + MRPapproveRRPapprovedummy, data = right_only)

# robust SEs
robust_delete_tweet_nocontrols_fullsample <- summ(delete_tweet_nocontrols_fullsample, robust = "HC1")
robust_delete_tweet_nocontrols_leftonly <- summ(delete_tweet_nocontrols_leftonly, robust = "HC1")
robust_delete_tweet_nocontrols_rightonly <- summ(delete_tweet_nocontrols_rightonly, robust = "HC1")

# convert regression outcomes to df for coefficient plot
robust_delete_tweet_nocontrols_fullsample_df <- broom::tidy(robust_delete_tweet_nocontrols_fullsample) %>%
  mutate(sample = "Full Sample") %>%
  mutate(measure = "Sanctioning") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_delete_tweet_nocontrols_fullsample_df)[1] <- "model"
names(robust_delete_tweet_nocontrols_fullsample_df)[7] <- "term"

robust_delete_tweet_nocontrols_leftonly_df <- broom::tidy(robust_delete_tweet_nocontrols_leftonly) %>%
  mutate(sample = "Left Wing Only") %>%
  mutate(measure = "Sanctioning") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_delete_tweet_nocontrols_leftonly_df)[1] <- "model"
names(robust_delete_tweet_nocontrols_leftonly_df)[7] <- "term"


robust_delete_tweet_nocontrols_rightonly_df <- broom::tidy(robust_delete_tweet_nocontrols_rightonly) %>%
  mutate(sample = "Right Wing Only") %>%
  mutate(measure = "Sanctioning") %>%
  filter(term != "(Intercept)")

# rename for the dwplot
names(robust_delete_tweet_nocontrols_rightonly_df)[1] <- "model"
names(robust_delete_tweet_nocontrols_rightonly_df)[7] <- "term"


# final joining of fig 1 - comparison against control in full sample
joined_models_nocontrols_fullsample <- rbind(robust_sign_petition_self_nocontrols_fullsample_df, robust_sensitive_item_nocontrols_fullsample_df, robust_petition_appropriate_self_nocontrols_fullsample_df, robust_sign_petition_others_nocontrols_fullsample_df, robust_petition_appropriate_others_nocontrols_fullsample_df, robust_delete_tweet_nocontrols_fullsample_df)

# Change a few labels
joined_models_nocontrols_fullsample$term[joined_models_nocontrols_fullsample$term == "Agreement with Sensitive Item"] <- "Agree Sensitive Item"
joined_models_nocontrols_fullsample$term[joined_models_nocontrols_fullsample$term == "Personal Willingness to Sign Petition"] <- "Sign Petition"
joined_models_nocontrols_fullsample$term[joined_models_nocontrols_fullsample$term == "Personal Views of Appropriateness of Signing"] <- "Perceived Appropriateness"


# Plot with the whole sample
fig_1 <- joined_models_nocontrols_fullsample %>%
  mutate(persuassion = ifelse(term %in% c(
    "Agree Sensitive Item",
    # "Personal Views of\nAppropriateness of\nSigning",
    "Sign Petition"
  ),
  "Preferences", "Norms Items"
  )) %>%
  mutate(facet_order = ifelse(persuassion == "Preferences", 0, 1)) %>%
  mutate(persuassion = fct_reorder(persuassion, facet_order)) %>%
  mutate(color_variable = ifelse(model %in% c("MRPapprovedummy", "MRPapproveRRPapprovedummy"),
    "Mainstream right approves", "Mainstream right does not approve"
  )) %>%
  mutate(sample_order = ifelse(model %in% c("MRPapprovedummy"), 3,
    ifelse(model %in% c("MRPapproveRRPapprovedummy"), 2,
      ifelse(model %in% c("MRPdisapproveRRPapprovedummy"), 1,
        0
      )
    )
  )) %>%
  mutate(model = fct_reorder(model, sample_order)) %>%
  mutate(term_order = ifelse(term %in% c("Sanctioning"), 0,
    ifelse(term %in% c("Normative Expectations"), 1,
      ifelse(term %in% c("Empirical Expectations"), 2,
        3
      )
    )
  )) %>%
  mutate(term = fct_reorder(term, term_order)) %>%
  mutate(ci_upper = estimate + 1.96 * std.error) %>%
  mutate(ci_lower = estimate - 1.96 * std.error) %>%
  mutate(ci_lower_90 = estimate - 1.645 * std.error) %>%
  mutate(ci_upper_90 = estimate + 1.645 * std.error) %>%
  ggplot(aes(term, estimate, shape = model, color = color_variable)) +
  geom_hline(aes(yintercept = 0), linetype = 2) +
  coord_flip(ylim = c(-0.5, 0.5)) +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), position = position_dodge(width = 0.75), size = 0.15) +
  geom_linerange(aes(ymin = ci_lower_90, ymax = ci_upper_90), position = position_dodge(width = 0.75), size = 0.5) +
  geom_point(size = 1, position = position_dodge(width = 0.75)) +
  scale_color_manual(values = c("black", "gray70")) +
  facet_grid(
    rows = vars(persuassion),
    scales = "free_y"
  ) +
  theme_minimal(base_size = 9) +
  theme(
    legend.title = element_blank(),
    panel.background = element_rect(fill = NA, color = "black")
  ) +
  scale_shape_discrete(labels = c(
    "RRP Approve",
    "MRP Disapprove and RRP Approve",
    "MRP Approve and RRP Approve",
    "MRP Approve"
  )) +
  guides(color = FALSE) +
  guides(shape = guide_legend(reverse = TRUE)) +
  ylab(" ") +
  xlab("")
# fig_1

ggsave("plots/fig_1.png", plot = fig_1, width = 17.5, height = 12, units = "cm")


# final joining of fig 2 - comparison against control in left wing, right wing and rrp samples
joined_models_nocontrols_partysamples <- rbind(
  robust_sensitive_item_nocontrols_leftonly_df, robust_sensitive_item_nocontrols_rightonly_df,
  robust_sign_petition_self_nocontrols_leftonly_df, robust_sign_petition_self_nocontrols_rightonly_df,
  robust_petition_appropriate_self_nocontrols_leftonly_df, robust_petition_appropriate_self_nocontrols_rightonly_df,
  robust_sign_petition_others_nocontrols_leftonly_df, robust_sign_petition_others_nocontrols_rightonly_df,
  robust_petition_appropriate_others_nocontrols_leftonly_df, robust_petition_appropriate_others_nocontrols_rightonly_df,
  robust_delete_tweet_nocontrols_leftonly_df, robust_delete_tweet_nocontrols_rightonly_df
)

# reorder to specify sequence in facet_wrap
joined_models_nocontrols_partysamples$sample <- factor(joined_models_nocontrols_partysamples$sample, # Reordering group factor levels
  levels = c("Right Wing Only", "Left Wing Only", "Radical Right Only")
)


# Change a few labels
joined_models_nocontrols_partysamples$term[joined_models_nocontrols_partysamples$term == "Agreement with Sensitive Item"] <- "Agree Sensitive Item"
joined_models_nocontrols_partysamples$term[joined_models_nocontrols_partysamples$term == "Personal Willingness to Sign Petition"] <- "Sign Petition"
joined_models_nocontrols_partysamples$term[joined_models_nocontrols_partysamples$term == "Personal Views of Appropriateness of Signing"] <- "Perceived Appropriateness"


# Final plot split by sample: Fig 2
fig_2 <- joined_models_nocontrols_partysamples %>%
  mutate(persuassion = ifelse(term %in% c(
    "Agree Sensitive Item",
    # "Personal Views of\nAppropriateness of\nSigning",
    "Sign Petition"
  ),
  "Preferences", "Norms Items"
  )) %>%
  mutate(facet_order = ifelse(persuassion == "Preferences", 0, 1)) %>%
  mutate(persuassion = fct_reorder(persuassion, facet_order)) %>%
  mutate(color_variable = ifelse(model %in% c("MRPapprovedummy", "MRPapproveRRPapprovedummy"),
    "Mainstream right approves", "Mainstream right does not approve"
  )) %>%
  mutate(sample_order = ifelse(model %in% c("MRPapprovedummy"), 3,
    ifelse(model %in% c("MRPapproveRRPapprovedummy"), 2,
      ifelse(model %in% c("MRPdisapproveRRPapprovedummy"), 1,
        0
      )
    )
  )) %>%
  mutate(model = fct_reorder(model, sample_order)) %>%
  mutate(term_order = ifelse(model %in% c("Sanctioning"), 3,
    ifelse(model %in% c("MRPapproveRRPapprovedummy"), 2,
      ifelse(model %in% c("MRPdisapproveRRPapprovedummy"), 1,
        0
      )
    )
  )) %>%
  mutate(term = fct_reorder(term, term_order)) %>%
  mutate(term_order = ifelse(term %in% c("Sanctioning"), 0,
    ifelse(term %in% c("Normative Expectations"), 1,
      ifelse(term %in% c("Empirical Expectations"), 2,
        3
      )
    )
  )) %>%
  mutate(term = fct_reorder(term, term_order)) %>%
  mutate(ci_upper = estimate + 1.96 * std.error) %>%
  mutate(ci_lower = estimate - 1.96 * std.error) %>%
  mutate(ci_lower_90 = estimate - 1.645 * std.error) %>%
  mutate(ci_upper_90 = estimate + 1.645 * std.error) %>%
  ggplot(aes(term, estimate, shape = model, color = color_variable)) +
  geom_hline(aes(yintercept = 0), linetype = 2) +
  coord_flip(ylim = c(-0.5, 0.5)) +
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), position = position_dodge(width = 0.75), size = 0.15) +
  geom_linerange(aes(ymin = ci_lower_90, ymax = ci_upper_90), position = position_dodge(width = 0.75), size = 0.5) +
  geom_point(size = 1, position = position_dodge(width = 0.75)) +
  scale_color_manual(values = c("black", "gray70")) +
  facet_grid(
    rows = vars(persuassion), cols = vars(sample),
    scales = "free_y"
  ) +
  theme_minimal(base_size = 9) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    panel.background = element_rect(fill = NA, color = "black")
  ) +
  scale_shape_discrete(labels = c(
    "RRP Approve",
    "MRP Disapprove and RRP Approve",
    "MRP Approve and RRP Approve",
    "MRP Approve"
  )) +
  guides(color = FALSE) +
  guides(shape = guide_legend(reverse = TRUE)) +
  ylab(" ") +
  xlab("")
# fig_2

ggsave("plots/fig_2.png", plot = fig_2, width = 20, height = 15, units = "cm")

# Check size of coefficients
# View(joined_models_nocontrols_partysamples)
# View(joined_models_nocontrols_partysamples[joined_models_nocontrols_partysamples$sample == "Left Wing Only"])
min(abs(joined_models_nocontrols_partysamples$estimate[joined_models_nocontrols_partysamples$sample == "Left Wing Only"]))
table(joined_models_nocontrols_partysamples$model[joined_models_nocontrols_partysamples$estimate == 0.00760726])

median(abs(joined_models_nocontrols_partysamples$estimate[joined_models_nocontrols_partysamples$sample == "Left Wing Only"]))
max(joined_models_nocontrols_partysamples$estimate[joined_models_nocontrols_partysamples$sample == "Left Wing Only"])

min(abs(joined_models_nocontrols_partysamples$estimate[joined_models_nocontrols_partysamples$sample == "Right Wing Only"]))
median(abs(joined_models_nocontrols_partysamples$estimate[joined_models_nocontrols_partysamples$sample == "Right Wing Only"]))
max(joined_models_nocontrols_partysamples$estimate[joined_models_nocontrols_partysamples$sample == "Right Wing Only"])

# Adjusting for multiple comparisons (for full sample analyses of Fig 1): Table E.1
joined_models_nocontrols_fullsample <- joined_models_nocontrols_fullsample %>%
  group_by(model) %>%
  mutate(
    p_bonferroni = round(p.adjust(p.value, method = "bonferroni", n = 6), 4),
    p_holm = round(p.adjust(p.value, method = "holm", n = 6), 4),
    p_BH = round(p.adjust(p.value, method = "BH", n = 6), 4)
  ) %>%
  ungroup()

adj_comps <- joined_models_nocontrols_fullsample %>%
  select(
    term, model, estimate, p.value,
    p_bonferroni, p_holm, p_BH
  )

adj_comps_tab <- xtable(adj_comps, digits = 3)

names(adj_comps_tab) <- c(
  "Outcome", "Treatment", "Coefficient", "P-value",
  "Corrected p-value (Bonferroni)",
  "Corrected p-value (Holm)",
  "Corrected p-value (Benjamini & Hochberg)"
)
print(adj_comps_tab, type = "html", file = "tables/table_e1.html")

# Multiple comparisons adjustments (for subsample analyses of fig 2): Tables E3 and E4
padj <- joined_models_nocontrols_partysamples %>%
  group_by(model, sample) %>%
  mutate(
    p_bonferroni = p.adjust(p.value, method = "bonferroni", n = 6),
    p_holm = p.adjust(p.value, method = "holm", n = 6),
    p_BH = p.adjust(p.value, method = "BH", n = 6)
  ) %>%
  ungroup() %>%
  group_by(sample)

adj_comps <- padj %>%
  select(
    sample, term, model, estimate, p.value,
    p_bonferroni, p_holm, p_BH
  ) %>%
  arrange(sample)

adj_comps_tab <- xtable(adj_comps, digits = 3)
names(adj_comps_tab) <- c(
  "Sample", "Outcome", "Treatment", "Coefficient",
  "P-value", "Corrected p-value (Bonferroni)",
  "Corrected p-value (Holm)",
  "Corrected p-value (Benjamini & Hochberg)"
)

print(adj_comps_tab, type = "html", file = "tables/tables_e3_e4.html")
